function [y,fy,d] = newton (f, df, y, a, b_f, b_df)

% Newton's Method
% ~~~~~~~~~~~~~~~
% INPUT:
%    f   => function name
%    df  => jacobian-function name
%    y   => initial guess
%    The following parameters are designed to fit the needs of my ODE solvers
%    a,b_f,b_df => (optional)
%                  a    => parameter for f,df (the same for both)
%                  b_f  => parameter for f
%                  b_df => parameter for df
%
% OUTPUT:
%    y   => zero of f
%    The following parameter is designed to fit the needs of my ODE solvers
%    fy  => if nargin=6, returns f(y) (close to zero, hopefully)
%           otherwise, returns the second output argument of f

if nargin == 3
   fy = feval(f,y);
   if norm(fy) > 1e-25
      d  = -feval(df,y) \ fy;
      while (norm(d)/(norm(y)+1e-5)) > 1e-8
         y  = y + d;
         fy = feval(f,y);
         d  = -feval(df,y) \ fy;
      end
   end
elseif nargin == 6
   [Fy,fy] = feval(f,y,b_f,a);
   if norm(Fy) > 1e-25
      d  = -feval(df,y,b_df,a) \ Fy;
      while (norm(d)/(norm(y)+1e-5)) > 1e-8
         y = y + d;
         [Fy,fy] = feval(f,y,b_f,a);
         d  = -feval(df,y,b_df,a) \ Fy;
      end
   end
else
   error('Wrong number of input arguments in function newton!');
end
